Notebook for cell type identity correlations using pairwise genes. Subsetting will be needed since for some of the species we only used the first cell type name (which is what we’re more interested in).
Setup chunk
Setup reticulate
knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")
[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/pallium_evo"
Load libraries
library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)
[1] FALSE
use_python(Sys.which("python"))
py_config()
python: /home/tpires/bin/miniconda3/bin/python
libpython: /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome: /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version: 3.8.3 (default, May 19 2020, 18:47:26) [GCC 7.3.0]
numpy: /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version: 1.18.5
NOTE: Python version was forced by RETICULATE_PYTHON
Load datasets
library(Seurat)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching SeuratObject
library(simspec)
library(harmony)
Loading required package: Rcpp
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
Registered S3 method overwritten by 'cli':
method from
print.boxx spatstat.geom
── Attaching packages ─────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.3 ✓ purrr 0.3.4
✓ tibble 3.1.2 ✓ stringr 1.4.0
✓ tidyr 1.1.3 ✓ forcats 0.5.0
✓ readr 1.4.0
── Conflicts ────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(ggplot2)
library(uwot)
Loading required package: Matrix
Attaching package: ‘Matrix’
The following objects are masked from ‘package:tidyr’:
expand, pack, unpack
Load metadatas
lizard_all = readRDS("data/expression/lizard_all_v3.RDS")
turtle_all = readRDS("data/expression/turtle_all_v3.RDS")
axolotl_all = readRDS(file = "data/expression/axolotl/palliumres07_harmony.RDS")
Add metadata
lizard_meta = read.csv("data/annotations/lizard_all_umeta.csv", header = T, row.names = 1)
turtle_meta = read.csv("data/annotations/turtle_all_umeta.csv", header = T, row.names = 1)
axolotl_meta = read.csv("data/annotations/axolotl_all_umeta.csv", header = T, row.names = 1)
Filter data (doublets and odd cells)
lizard_all = AddMetaData(lizard_all, metadata = lizard_meta)
turtle_all = AddMetaData(turtle_all, metadata = turtle_meta)
axolotl_all = AddMetaData(axolotl_all, metadata = axolotl_meta)
Put Seurats in list
turtle_all = turtle_all[,turtle_all$classes!="doublets"]
axolotl_all = axolotl_all[,!axolotl_all$cellclusters %in% c("GABAergic neuron", "ETV1+ neuron", "neuron")]
Load necessary orthology data
srat_list = list("paintedturtle" = turtle_all, "axolotl" = axolotl_all, "lizard" = lizard_all)
Function for ortholog matching
ort_alt_ens_list = readRDS(file = "data/eggNOG/ort_alt_ens_list.RDS")
ort_alt_egg_list = readRDS(file = "data/eggNOG/ort_alt_egg_list.RDS")
ort_list_use = c(ort_alt_ens_list,
ort_alt_egg_list[which(grepl("axolotl", names(ort_alt_egg_list)) |
grepl("lizard", names(ort_alt_egg_list)))])
bfmats = ort_alt_egg_list[which(grepl("zebrafinch", names(ort_alt_egg_list)))]
names(bfmats) = gsub("zebrafinch", "bengalesefinch", names(bfmats))
ort_list_use = c(ort_list_use, bfmats)
Function for subsetting Seurats
# get orthologs for all species pairs
ortMatch = function(ort_l, sp_vec, ort_scope = "ortholog_one2one"){
sp_p_mat = combn(sp_vec, 2)
for(i in 1:ncol(sp_p_mat)){
sp_pair = sp_p_mat[,i]
ort_tab = ort_l[[paste0(sp_pair, collapse = ".vs.")]]
if(!is.null(ort_scope)){
ort_tab = ort_tab[ort_tab$homology.type==ort_scope,]
}
colnames(ort_tab) = gsub("genename", "", colnames(ort_tab))
if(i==1){
ort_all = unique(ort_tab[,c(2,4)])
} else{
ort_all = unique(merge(ort_all, unique(ort_tab[,c(2,4)]),
by = intersect(colnames(ort_all), colnames(ort_tab))))
}
}
return(ort_all)
}
Get ortholog-matched Seurats
seuratOrthologs = function(s_l, o_all){
# only orthologs present in the Seurat objects
for(sp in names(s_l)){
o_all = o_all[o_all[,sp] %in% rownames(s_l[[sp]]),]
}
# get a unique name for the genes
allg = unlist(o_all)
jointnames = c()
rep_n = list()
for(i in 1:nrow(o_all)){
nname = unlist(lapply(o_all[i,], function(x) sum(allg==x)))
g = unlist(o_all[i,names(nname[nname==1])])
# if there are acceptable unique names
if(any(nname==1) & (!is.null(g) & !all(grepl("..", g, fixed = T)))){ # no weird axololt annot
g = g[!grepl("..", g, fixed = T)] # no weird axololt annot
g = c(g[!grepl("LOC", g)], g[grepl("LOC", g)])[1] #prioritize non-LOC
jointnames = c(jointnames, g)
rep_n[[g]] = 1
} else{ # we'll add a number to those names
g = unlist(o_all[i,])
g = g[!grepl("..", g, fixed = T)] # no weird axololt annot
g = c(g[!grepl("LOC", g)], g[grepl("LOC", g)])[1] #prioritize non-LOC
if(g %in% names(rep_n)){
gn = rep_n[[g]]
rep_n[[g]] = gn+1
jointnames = c(jointnames, paste0(g, "-", gn+1))
} else{
rep_n[[g]] = 1
jointnames = c(jointnames, paste0(g, "-", 1))
}
}
}
# filter the Seurats by the orthologs
## this requires redoing the Seurat with new names, since some genes will be duplicated
for(sp in names(s_l)){
cnts = s_l[[sp]]@assays$RNA@counts[o_all[,sp],]
rownames(cnts) = jointnames
m = s_l[[sp]]@meta.data
s_l[[sp]] = CreateSeuratObject(cnts, project = sp, meta.data = m)
}
return(s_l)
}
Normalise data
for(sp in names(srat_ort)){
print(sp)
srat_ort[[sp]] = suppressWarnings(SCTransform(srat_ort[[sp]], do.correct.umi = T, verbose = F,
seed.use = 1, vars.to.regress = "nCount_RNA",
variable.features.rv.th = 1, return.only.var.genes = F,
variable.features.n = NULL))
srat_ort1[[sp]] = suppressWarnings(SCTransform(srat_ort1[[sp]], do.correct.umi = T, verbose = F,
seed.use = 1, vars.to.regress = "nCount_RNA",
variable.features.rv.th = 1, return.only.var.genes = F,
variable.features.n = NULL))
}
Seurats with all merged
for(sp in names(srat_ort)){
print(sp)
srat_ort[[sp]] = suppressWarnings(SCTransform(srat_ort[[sp]], do.correct.umi = T, verbose = F,
seed.use = 1, vars.to.regress = "nCount_RNA",
variable.features.rv.th = 1, return.only.var.genes = F,
variable.features.n = NULL))
srat_ort1[[sp]] = suppressWarnings(SCTransform(srat_ort1[[sp]], do.correct.umi = T, verbose = F,
seed.use = 1, vars.to.regress = "nCount_RNA",
variable.features.rv.th = 1, return.only.var.genes = F,
variable.features.n = NULL))
}
[1] "paintedturtle"
[1] "axolotl"
[1] "lizard"
Integration with all ortholog
srat_ort_all = Reduce(merge, srat_ort)
srat_ort_all = AddMetaData(srat_ort_all,
metadata = unlist(lapply(names(srat_ort),
function(x) rep(x, ncol(srat_ort[[x]])))),
col.name = "species")
srat_ort1_all = Reduce(merge, srat_ort1)
srat_ort1_all = AddMetaData(srat_ort1_all,
metadata = unlist(lapply(names(srat_ort1),
function(x) rep(x, ncol(srat_ort1[[x]])))),
col.name = "species")
Integration with one2one ortholog
feat_list = list("all" = rownames(srat_ort1_all@assays$SCT@scale.data),
"allvar" = unique(unlist(lapply(srat_ort1, VariableFeatures))),
"commonvar" = Reduce(intersect, lapply(srat_ort1, VariableFeatures)))
varindat = lapply(srat_ort1, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort1)]
srat_rpca_one_l = list()
for(n in names(feat_list)){
print(n)
integr_feat = feat_list[[n]]
# calculate possible missing Pearson residuals for SCTransform
srat_ort1_int = PrepSCTIntegration(srat_ort1, anchor.features = integr_feat, verbose = T)
srat_ort1_int <- lapply(X = srat_ort1_int, FUN = function(x) {
x <- ScaleData(x, features = integr_feat, verbose = FALSE)
x <- RunPCA(x, features = integr_feat, verbose = FALSE)
})
# finding the anchors for integration
all_cell_anchors = FindIntegrationAnchors(srat_ort1_int, normalization.method = "SCT",
anchor.features = integr_feat, reduction = "rpca",
assay = rep("SCT", length(srat_ort1_int)),
dims = 1:50, verbose = T)
# actual data integration
allgenes = rownames(srat_ort1_int[[1]]@assays$RNA@counts)
all_cell_integr = IntegrateData(all_cell_anchors, normalization.method = "SCT", dims = 1:50,
verbose = T, features.to.integrate = allgenes)
# run PCA and UMAP to see how it looks
all_cell_integr = RunPCA(all_cell_integr, verbose = F)
all_cell_integr = RunUMAP(all_cell_integr, dims = 1:30)
all_cell_integr = AddMetaData(all_cell_integr,
metadata = unlist(lapply(names(srat_ort1),
function(x) rep(x, ncol(srat_ort1[[x]])))),
col.name = "species")
srat_rpca_one_l[[n]] = all_cell_integr
}
Integration with all ortholog
feat_list = list("all" = rownames(srat_ort_all@assays$SCT@scale.data),
"allvar" = unique(unlist(lapply(srat_ort, VariableFeatures))),
"commonvar" = Reduce(intersect, lapply(srat_ort, VariableFeatures)))
varindat = lapply(srat_ort, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort)]
srat_harm_all_l = list()
for(n in names(feat_list)){
print(n)
# Renormalise
integr_feat = feat_list[[n]]
VariableFeatures(srat_ort_all) = integr_feat
srat_ort_all_int = ScaleData(srat_ort_all, features = integr_feat, use.umi = T,
do.scale = F, verbose = FALSE)
srat_ort_all_int = RunPCA(srat_ort_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
features = integr_feat)
# Run Harmony
srat_ort_all_int = RunHarmony(srat_ort_all_int, "species", tau = 30,
plot_convergence = F, assay.use = "SCT")
# Run UMAP on Harmony dimensions
srat_ort_all_int = RunUMAP(srat_ort_all_int, reduction = "harmony", dims = 1:10)
srat_harm_all_l[[n]] = srat_ort_all_int
}
Integration with one2one ortholog
feat_list = list("all" = rownames(srat_ort1_all@assays$SCT@scale.data),
"allvar" = unique(unlist(lapply(srat_ort1, VariableFeatures))),
"commonvar" = Reduce(intersect, lapply(srat_ort1, VariableFeatures)))
varindat = lapply(srat_ort1, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort1)]
srat_harm_one_l = list()
for(n in names(feat_list)){
print(n)
# Renormalise
integr_feat = feat_list[[n]]
VariableFeatures(srat_ort1_all) = integr_feat
srat_ort1_all_int = ScaleData(srat_ort1_all, features = integr_feat, use.umi = T,
do.scale = F, verbose = FALSE)
srat_ort1_all_int = RunPCA(srat_ort1_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
features = integr_feat)
# Run Harmony
srat_ort1_all_int = RunHarmony(srat_ort1_all_int, "species", tau = 30,
plot_convergence = F, assay.use = "SCT")
# Run UMAP on Harmony dimensions
srat_ort1_all_int = RunUMAP(srat_ort1_all_int, reduction = "harmony", dims = 1:10)
srat_harm_one_l[[n]] = srat_ort1_all_int
}
Integration with all ortholog
feat_list = list("all" = rownames(srat_ort_all@assays$SCT@scale.data),
"allvar" = unique(unlist(lapply(srat_ort, VariableFeatures))),
"commonvar" = Reduce(intersect, lapply(srat_ort, VariableFeatures)))
varindat = lapply(srat_ort, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort)]
srat_css_all_l = list()
for(n in names(feat_list)){
print(n)
integr_feat = feat_list[[n]]
VariableFeatures(srat_ort_all) = integr_feat
srat_ort_all_int = ScaleData(srat_ort_all, features = integr_feat, use.umi = T,
do.scale = F, verbose = FALSE)
srat_ort_all_int = RunPCA(srat_ort_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
features = integr_feat)
# Run CSS
srat_ort_all_int = cluster_sim_spectrum(srat_ort_all_int, label_tag="species", redo_pca = T,
dims_use = 1:30, var_genes = integr_feat)
# Run UMAP on CSS dimensions
srat_ort_all_int = RunUMAP(srat_ort_all_int, reduction="css", reduction.name="umap_css",
dims = 1:ncol(Embeddings(srat_ort_all_int,"css")), reduction.key="UMAPCSS_")
srat_ort_all_int = SetIdent(srat_ort_all_int, value = "species")
srat_css_all_l[[n]] = srat_ort_all_int
}
Integration with one2one ortholog
feat_list = list("all" = rownames(srat_ort_all@assays$SCT@scale.data),
"allvar" = unique(unlist(lapply(srat_ort, VariableFeatures))),
"commonvar" = Reduce(intersect, lapply(srat_ort, VariableFeatures)))
varindat = lapply(srat_ort, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort)]
srat_css_all_l = list()
for(n in names(feat_list)){
print(n)
integr_feat = feat_list[[n]]
VariableFeatures(srat_ort_all) = integr_feat
srat_ort_all_int = ScaleData(srat_ort_all, features = integr_feat, use.umi = T,
do.scale = F, verbose = FALSE)
srat_ort_all_int = RunPCA(srat_ort_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
features = integr_feat)
# Run CSS
srat_ort_all_int = cluster_sim_spectrum(srat_ort_all_int, label_tag="species", redo_pca = T,
dims_use = 1:30, var_genes = integr_feat)
# Run UMAP on CSS dimensions
srat_ort_all_int = RunUMAP(srat_ort_all_int, reduction="css", reduction.name="umap_css",
dims = 1:ncol(Embeddings(srat_ort_all_int,"css")), reduction.key="UMAPCSS_")
srat_ort_all_int = SetIdent(srat_ort_all_int, value = "species")
srat_css_all_l[[n]] = srat_ort_all_int
}
[1] "all"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 16415
Number of edges: 611282
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9237
Number of communities: 20
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4187
Number of edges: 160419
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8918
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 18664
Number of edges: 619386
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8942
Number of communities: 15
Elapsed time: 2 seconds
1 singletons identified. 14 final clusters.
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.
23:28:48 UMAP embedding parameters a = 0.9922 b = 1.112
23:28:48 Read 39266 rows and found 45 numeric columns
23:28:48 Using Annoy for neighbor search, n_neighbors = 30
23:28:48 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:28:57 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd14969e74
23:28:57 Searching Annoy index using 1 thread, search_k = 3000
23:29:10 Annoy recall = 100%
23:29:13 Commencing smooth kNN distance calibration using 1 thread
23:29:17 Initializing from normalized Laplacian + noise
23:29:21 Commencing optimization for 200 epochs, with 1762176 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:30:11 Optimization finished
[1] "allvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 16415
Number of edges: 613106
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9233
Number of communities: 19
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4187
Number of edges: 160511
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8930
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 18664
Number of edges: 615768
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8890
Number of communities: 14
Elapsed time: 2 seconds
1 singletons identified. 13 final clusters.
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.
23:40:23 UMAP embedding parameters a = 0.9922 b = 1.112
23:40:23 Read 39266 rows and found 43 numeric columns
23:40:23 Using Annoy for neighbor search, n_neighbors = 30
23:40:23 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:40:29 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd16b007fa
23:40:29 Searching Annoy index using 1 thread, search_k = 3000
23:40:42 Annoy recall = 100%
23:40:44 Commencing smooth kNN distance calibration using 1 thread
23:40:48 Initializing from normalized Laplacian + noise
23:40:52 Commencing optimization for 200 epochs, with 1764938 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:41:42 Optimization finished
[1] "commonvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 16415
Number of edges: 612295
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9173
Number of communities: 16
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4187
Number of edges: 169158
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8871
Number of communities: 10
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 18664
Number of edges: 642681
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8838
Number of communities: 14
Elapsed time: 2 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.
23:44:05 UMAP embedding parameters a = 0.9922 b = 1.112
23:44:05 Read 39266 rows and found 40 numeric columns
23:44:05 Using Annoy for neighbor search, n_neighbors = 30
23:44:05 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:44:11 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd148aac04
23:44:11 Searching Annoy index using 1 thread, search_k = 3000
23:44:24 Annoy recall = 100%
23:44:25 Commencing smooth kNN distance calibration using 1 thread
23:44:29 Initializing from normalized Laplacian + noise
23:44:33 Commencing optimization for 200 epochs, with 1796060 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:45:23 Optimization finished
rPCA + Harmony, all ortholog
feat_list = list("all" = rownames(srat_ort1_all@assays$SCT@scale.data),
"allvar" = unique(unlist(lapply(srat_ort1, VariableFeatures))),
"commonvar" = Reduce(intersect, lapply(srat_ort1, VariableFeatures)))
varindat = lapply(srat_ort1, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort1)]
srat_css_one_l = list()
for(n in names(feat_list)){
print(n)
integr_feat = feat_list[[n]]
VariableFeatures(srat_ort1_all) = integr_feat
srat_ort1_all_int = ScaleData(srat_ort1_all, features = integr_feat, use.umi = T,
do.scale = F, verbose = FALSE)
srat_ort1_all_int = RunPCA(srat_ort1_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
features = integr_feat)
# Run CSS
srat_ort1_all_int = cluster_sim_spectrum(srat_ort1_all_int, label_tag="species", redo_pca = T,
dims_use = 1:30, var_genes = integr_feat)
# Run UMAP on CSS dimensions
srat_ort1_all_int = RunUMAP(srat_ort1_all_int, reduction="css", reduction.name="umap_css",
dims = 1:ncol(Embeddings(srat_ort1_all_int,"css")),
reduction.key="UMAPCSS_")
srat_ort1_all_int = SetIdent(srat_ort1_all_int, value = "species")
srat_css_one_l[[n]] = srat_ort1_all_int
}
[1] "all"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 16415
Number of edges: 640464
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9261
Number of communities: 20
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4187
Number of edges: 163557
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8957
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 18664
Number of edges: 626082
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8882
Number of communities: 13
Elapsed time: 3 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.
23:55:34 UMAP embedding parameters a = 0.9922 b = 1.112
23:55:34 Read 39266 rows and found 44 numeric columns
23:55:34 Using Annoy for neighbor search, n_neighbors = 30
23:55:34 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:55:40 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd3ecb81da
23:55:40 Searching Annoy index using 1 thread, search_k = 3000
23:55:53 Annoy recall = 100%
23:55:55 Commencing smooth kNN distance calibration using 1 thread
23:55:59 Initializing from normalized Laplacian + noise
23:56:03 Commencing optimization for 200 epochs, with 1777116 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:56:53 Optimization finished
[1] "allvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 16415
Number of edges: 645110
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9265
Number of communities: 21
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4187
Number of edges: 165059
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8960
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 18664
Number of edges: 622418
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8864
Number of communities: 12
Elapsed time: 2 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.
00:04:34 UMAP embedding parameters a = 0.9922 b = 1.112
00:04:34 Read 39266 rows and found 44 numeric columns
00:04:34 Using Annoy for neighbor search, n_neighbors = 30
00:04:34 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:04:40 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd505d4a3f
00:04:40 Searching Annoy index using 1 thread, search_k = 3000
00:04:54 Annoy recall = 100%
00:04:55 Commencing smooth kNN distance calibration using 1 thread
00:04:59 Initializing from normalized Laplacian + noise
00:05:03 Commencing optimization for 200 epochs, with 1777102 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:05:53 Optimization finished
[1] "commonvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 16415
Number of edges: 623330
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.9174
Number of communities: 17
Elapsed time: 1 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 4187
Number of edges: 172270
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8892
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 18664
Number of edges: 638406
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8901
Number of communities: 13
Elapsed time: 2 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.
00:07:52 UMAP embedding parameters a = 0.9922 b = 1.112
00:07:52 Read 39266 rows and found 41 numeric columns
00:07:52 Using Annoy for neighbor search, n_neighbors = 30
00:07:52 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:07:58 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd60b3d33d
00:07:58 Searching Annoy index using 1 thread, search_k = 3000
00:08:13 Annoy recall = 100%
00:08:14 Commencing smooth kNN distance calibration using 1 thread
00:08:18 Initializing from normalized Laplacian + noise
00:08:21 Commencing optimization for 200 epochs, with 1810424 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:09:11 Optimization finished
rPCA + Harmony, one2one ortholog
srat_rpca_harm_one_l = list()
for(n in names(srat_rpca_one_l)){
print(n)
# Run Harmony
cell_integr = RunHarmony(srat_rpca_one_l[[n]], "species", tau = 30,
plot_convergence = F, assay.use = "SCT")
# Run UMAP on Harmony dimensions
cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
srat_rpca_harm_one_l[[n]] = cell_integr
}
CSS + Harmony, all ortholog
srat_css_harm_all_l = list()
for(n in names(srat_css_all_l)){
print(n)
# Run Harmony
cell_integr = RunHarmony(srat_css_all_l[[n]], "species", tau = 30, reduction = "css",
dims.use = 1:ncol(Embeddings(srat_css_all_l[[n]],"css")),
plot_convergence = F, assay.use = "SCT")
# Run UMAP on Harmony dimensions
cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
srat_css_harm_all_l[[n]] = cell_integr
}
CSS + Harmony, one2one ortholog
srat_css_harm_all_l = list()
for(n in names(srat_css_all_l)){
print(n)
# Run Harmony
cell_integr = RunHarmony(srat_css_all_l[[n]], "species", tau = 30, reduction = "css",
dims.use = 1:ncol(Embeddings(srat_css_all_l[[n]],"css")),
plot_convergence = F, assay.use = "SCT")
# Run UMAP on Harmony dimensions
cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
srat_css_harm_all_l[[n]] = cell_integr
}
[1] "all"
Harmony 1/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 7 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:11:22 UMAP embedding parameters a = 0.9922 b = 1.112
00:11:22 Read 39266 rows and found 10 numeric columns
00:11:22 Using Annoy for neighbor search, n_neighbors = 30
00:11:22 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:11:27 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd3ea2902b
00:11:27 Searching Annoy index using 1 thread, search_k = 3000
00:11:43 Annoy recall = 100%
00:11:45 Commencing smooth kNN distance calibration using 1 thread
00:11:48 Initializing from normalized Laplacian + noise
00:11:50 Commencing optimization for 200 epochs, with 1651614 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:12:39 Optimization finished
[1] "allvar"
Quick-TRANSfer stage steps exceeded maximum (= 1963300)Harmony 1/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 8/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 9/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 10/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:14:46 UMAP embedding parameters a = 0.9922 b = 1.112
00:14:46 Read 39266 rows and found 10 numeric columns
00:14:46 Using Annoy for neighbor search, n_neighbors = 30
00:14:46 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:14:51 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd41e214cd
00:14:51 Searching Annoy index using 1 thread, search_k = 3000
00:15:06 Annoy recall = 100%
00:15:07 Commencing smooth kNN distance calibration using 1 thread
00:15:10 Initializing from normalized Laplacian + noise
00:15:12 Commencing optimization for 200 epochs, with 1645504 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:16:01 Optimization finished
[1] "commonvar"
Harmony 1/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 7 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:17:37 UMAP embedding parameters a = 0.9922 b = 1.112
00:17:37 Read 39266 rows and found 10 numeric columns
00:17:37 Using Annoy for neighbor search, n_neighbors = 30
00:17:37 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:17:42 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd5e63e6fa
00:17:42 Searching Annoy index using 1 thread, search_k = 3000
00:17:58 Annoy recall = 100%
00:17:59 Commencing smooth kNN distance calibration using 1 thread
00:18:03 Initializing from normalized Laplacian + noise
00:18:05 Commencing optimization for 200 epochs, with 1670228 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:18:54 Optimization finished
Save
srat_css_harm_one_l = list()
for(n in names(srat_css_one_l)){
print(n)
# Run Harmony
cell_integr = RunHarmony(srat_css_one_l[[n]], "species", tau = 30, reduction = "css",
dims.use = 1:ncol(Embeddings(srat_css_one_l[[n]],"css")),
plot_convergence = F, assay.use = "SCT")
# Run UMAP on Harmony dimensions
cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
srat_css_harm_one_l[[n]] = cell_integr
}
[1] "all"
Harmony 1/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:20:34 UMAP embedding parameters a = 0.9922 b = 1.112
00:20:34 Read 39266 rows and found 10 numeric columns
00:20:34 Using Annoy for neighbor search, n_neighbors = 30
00:20:34 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:20:40 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd293b9c4c
00:20:40 Searching Annoy index using 1 thread, search_k = 3000
00:20:53 Annoy recall = 100%
00:20:55 Commencing smooth kNN distance calibration using 1 thread
00:20:58 Initializing from normalized Laplacian + noise
00:21:00 Commencing optimization for 200 epochs, with 1654888 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:21:49 Optimization finished
[1] "allvar"
Harmony 1/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 8/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 9/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 9 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:23:51 UMAP embedding parameters a = 0.9922 b = 1.112
00:23:51 Read 39266 rows and found 10 numeric columns
00:23:51 Using Annoy for neighbor search, n_neighbors = 30
00:23:51 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:23:56 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd4c0ca239
00:23:56 Searching Annoy index using 1 thread, search_k = 3000
00:24:10 Annoy recall = 100%
00:24:11 Commencing smooth kNN distance calibration using 1 thread
00:24:14 Initializing from normalized Laplacian + noise
00:24:16 Commencing optimization for 200 epochs, with 1651456 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:25:05 Optimization finished
[1] "commonvar"
did not converge in 25 iterationsHarmony 1/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:26:31 UMAP embedding parameters a = 0.9922 b = 1.112
00:26:31 Read 39266 rows and found 10 numeric columns
00:26:31 Using Annoy for neighbor search, n_neighbors = 30
00:26:31 Building Annoy index with metric = cosine, n_trees = 50
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:26:36 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd88888be
00:26:36 Searching Annoy index using 1 thread, search_k = 3000
00:26:51 Annoy recall = 100%
00:26:52 Commencing smooth kNN distance calibration using 1 thread
00:26:55 Initializing from normalized Laplacian + noise
00:26:57 Commencing optimization for 200 epochs, with 1660784 positive edges
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:27:46 Optimization finished
Plots for all genes
#saveRDS(srat_rpca_all_l, file = "./data/integration/srat_rpca_all_l.RDS")
#saveRDS(srat_rpca_one_l, file = "./data/integration/srat_rpca_one_l.RDS")
#saveRDS(srat_harm_all_l, file = "./data/integration/srat_harm_all_l.RDS")
#saveRDS(srat_harm_one_l, file = "./data/integration/srat_harm_one_l.RDS")
saveRDS(srat_css_all_l, file = "./data/integration/srat_css_all_l.RDS")
saveRDS(srat_css_one_l, file = "./data/integration/srat_css_one_l.RDS")
#saveRDS(srat_rpca_harm_all_l, file = "./data/integration/srat_rpca_harm_all_l.RDS")
#saveRDS(srat_rpca_harm_one_l, file = "./data/integration/srat_rpca_harm_one_l.RDS")
saveRDS(srat_css_harm_all_l, file = "./data/integration/srat_css_harm_all_l.RDS")
saveRDS(srat_css_harm_one_l, file = "./data/integration/srat_css_harm_one_l.RDS")
Plots for one2one genes
for(n in names(srat_rpca_all_l)){
plt1 = DimPlot(srat_rpca_all_l[[n]], reduction = "umap", group.by = "species")+
labs(title = paste0("rPCA, all orthologs, ", n))+
DimPlot(srat_rpca_all_l[[n]], reduction = "umap", group.by = "classes")+
labs(title = paste0("rPCA, all orthologs, ", n))
print(plt1)
plt2 = DimPlot(srat_css_all_l[[n]], reduction = "umap_css", group.by = "species")+
labs(title = paste0("CSS, all orthologs, ", n))+
DimPlot(srat_css_all_l[[n]], reduction = "umap_css", group.by = "classes")+
labs(title = paste0("CSS, all orthologs, ", n))
print(plt2)
plt3 = DimPlot(srat_harm_all_l[[n]], reduction = "umap", group.by = "species")+
labs(title = paste0("Harmony, all orthologs, ", n))+
DimPlot(srat_harm_all_l[[n]], reduction = "umap", group.by = "classes")+
labs(title = paste0("Harmony, all orthologs, ", n))
print(plt3)
plt4 = DimPlot(srat_rpca_harm_all_l[[n]], reduction = "umap", group.by = "species")+
labs(title = paste0("rPCA+Harmony, all orthologs, ", n))+
DimPlot(srat_rpca_harm_all_l[[n]], reduction = "umap", group.by = "classes")+
labs(title = paste0("rPCA+Harmony, all orthologs, ", n))
print(plt4)
plt5 = DimPlot(srat_css_harm_all_l[[n]], reduction = "umap", group.by = "species")+
labs(title = paste0("CSS+Harmony, all orthologs, ", n))+
DimPlot(srat_css_harm_all_l[[n]], reduction = "umap", group.by = "classes")+
labs(title = paste0("CSS+Harmony, all orthologs, ", n))
print(plt5)
}